#include<Matrix/LB_Sparse_Matrix.h>



void lb_sparse_matrix_free(LB_Sparse_Matrix*lsm)
{
    // printf("here\n");

    if(lsm==NULL){return;}
    // printf("here\n");

    for(int i=0;i<lsm->rows;i++)
    {
        SAFE_FREE(lsm->data[i]);
    }
    // printf("here\n");
    SAFE_FREE(lsm->data);
    SAFE_FREE(lsm->rowSizes);
    free(lsm);
}
void lb_sparse_matrix_resize(LB_Sparse_Matrix*lsm,int rows,int cols)
{

    SAFE_FREE(lsm->rowSizes);
    lsm->rowSizes=(int*)malloc(sizeof(int)*rows);
    memset(lsm->rowSizes,0,sizeof(int)*rows);
    if(lsm->data!=NULL)
    {
        for(int i=0;i<lsm->rows;i++)
        {
            SAFE_FREE(lsm->data[i]);
        }
        free(lsm->data);
        lsm->data=NULL;
    }

    lsm->data=(SparseMatrixEntry**)malloc(sizeof(SparseMatrixEntry*)*rows);
    memset(lsm->data,0,sizeof(SparseMatrixEntry*)*rows);
    lsm->rows=rows;lsm->cols=cols;
    lsm->prop=NULL;
}

void lb_sparse_matrix_copy_data(LB_Sparse_Matrix* lsm1,
    LB_Sparse_Matrix* lsm2)
{
    if(lsm1==NULL||lsm2==NULL){return;}
    lb_sparse_matrix_resize(lsm1,lsm2->rows,lsm2->cols);
    for(int i=0;i<lsm1->rows;i++)
    {
        lb_sparse_matrix_setrowsize(lsm1,i,lsm2->rowSizes[i]);
        memmove(lsm1->data[i],lsm2->data[i],
            sizeof(double)*lsm2->rowSizes[i]);
    }
}
LB_Sparse_Matrix* lb_sparse_matrix_new_copy(LB_Sparse_Matrix* lsm)
{

    LB_Sparse_Matrix* result = (LB_Sparse_Matrix*)malloc(sizeof(LB_Sparse_Matrix));
    lb_sparse_matrix_initn(result, lsm->rows, lsm->cols);
    for(int i=0;i<lsm->rows;i++)
    {
        lb_sparse_matrix_setrowsize(result,i,lsm->rowSizes[i]);
        memmove(result->data[i],lsm->data[i],
            sizeof(double)*lsm->rowSizes[i]);
    }

    return result;
}

void lb_sparse_matrix_multply_scala(LB_Sparse_Matrix* lsm,
    double fac,LB_Sparse_Matrix* re)
{
    if(lsm==NULL||re==NULL){return;}

    if(re->rows!=lsm->rows||re->cols!=lsm->cols)
    {
        // re->rows=lsm->rows;re->cols=lsm->cols;

        lb_sparse_matrix_resize(re,lsm->rows,lsm->cols);
    }

    for(int i=0;i<lsm->rows;i++)
    {
        if(re!=lsm)
        {
            re->rowSizes[i]=lsm->rowSizes[i];
            re->data[i]=(SparseMatrixEntry*)realloc(re->data[i],
                sizeof(SparseMatrixEntry)*re->rowSizes[i]);
        }

        for (int k = 0; k < lsm->rowSizes[i]; k++)
        {
            SparseMatrixEntry entry = lsm->data[i][k];
            SparseMatrixEntry entry1 = re->data[i][k];
            entry1.v=entry.v*fac;
        }
    }
}

void lb_sparse_matrix_multiply(LB_Sparse_Matrix*result,
    LB_Sparse_Matrix* lsm1,
    LB_Sparse_Matrix* lsm2)
{
    if (lsm1 == NULL || lsm2 == NULL || lsm1->cols != lsm2->rows) {
        return ; // 矩阵无效或维度不匹配
    }

    // 创建结果矩阵
    // LB_Sparse_Matrix* result = (LB_Sparse_Matrix*)malloc(sizeof(LB_Sparse_Matrix));
    // lb_sparse_matrix_initn(result, lsm1->rows, lsm2->cols);
    lb_sparse_matrix_resize(result, lsm1->rows, lsm2->cols);

    // int* nonZeroCounts =(int*)malloc(sizeof(int)*lsm1->rows);

    double* temp_value=(double*)malloc(sizeof(double)*lsm2->cols);

    for (int i = 0; i < lsm1->rows; i++) {
        memset(temp_value,0,sizeof(double)*lsm2->cols);
        for (int k = 0; k < lsm1->rowSizes[i]; k++)
        {
            SparseMatrixEntry entry = lsm1->data[i][k];
            for(int l=0;l<lsm2->rowSizes[entry.i];l++)
            {
                SparseMatrixEntry entry2 = lsm2->data[entry.i][l];
                temp_value[entry2.i]+=(entry.v*entry2.v);
            }
        }

        for(int j=0;j<lsm2->cols;j++)
        {
            if(fabs(temp_value[j])>1e-10)
            {
                result->rowSizes[i]++;
            }
        }
        result->data[i]=(SparseMatrixEntry*)malloc(sizeof(SparseMatrixEntry) *result->rowSizes[i]);
        int temp_sum=0;
        for(int j=0;j<lsm2->cols;j++)
        {
            if(fabs(temp_value[j])>1e-10)
            {
                result->data[i][temp_sum].i=j;
                result->data[i][temp_sum].v=temp_value[j];
                temp_sum++;
            }
        }
    }
    free(temp_value);
    // free(nonZeroCounts);
    // return result;
}


// result 可以等于lsm1 或者lsm2
//
void lb_sparse_matrix_plus(LB_Sparse_Matrix* result,
    LB_Sparse_Matrix* lsm1,
    LB_Sparse_Matrix* lsm2)
{
    if (lsm1 == NULL || lsm2 == NULL || lsm1->rows != lsm2->rows||
        lsm1->cols != lsm2->cols) {return ;}

    // 创建结果矩阵
    // LB_Sparse_Matrix* result = (LB_Sparse_Matrix*)malloc(sizeof(LB_Sparse_Matrix));
    // lb_sparse_matrix_initn(result, lsm1->rows, lsm2->cols);
    if(lsm1->rows!=result->rows||lsm1->cols!=result->cols)
    {
        lb_sparse_matrix_resize(result, lsm1->rows, lsm2->cols);
    }
    double* temp_value=(double*)malloc(sizeof(double)*lsm2->cols);

    for (int i = 0; i < lsm1->rows; i++)
    {
        memset(temp_value,0,sizeof(double)*lsm1->cols);
        // result->rowSizes[i]=0;
        for(int k=0;k<lsm1->rowSizes[i];k++)
        {
            SparseMatrixEntry entry = lsm1->data[i][k];
            temp_value[entry.i]+=entry.v;
        }

        for(int k=0;k<lsm2->rowSizes[i];k++)
        {
            SparseMatrixEntry entry = lsm2->data[i][k];
            temp_value[entry.i]+=entry.v;
        }
        result->rowSizes[i]=0;
        for(int j=0;j<lsm1->cols;j++)
        {
            if(fabs(temp_value[j])>1e-10)
            {
                result->rowSizes[i]++;
            }
        }
        SAFE_FREE(result->data[i]);
        result->data[i]=(SparseMatrixEntry*)malloc(
            sizeof(SparseMatrixEntry) *result->rowSizes[i]);
        int temp_sum=0;
        for(int j=0;j<lsm2->cols;j++)
        {
            if(fabs(temp_value[j])>1e-10)
            {
                result->data[i][temp_sum].i=j;
                result->data[i][temp_sum].v=temp_value[j];
                temp_sum++;
            }
        }
    }
    free(temp_value);
    // return result;
}

// result 可以等于lsm1 或者lsm2
//
void lb_sparse_matrix_minus(LB_Sparse_Matrix* result,
    LB_Sparse_Matrix* lsm1,
    LB_Sparse_Matrix* lsm2)
{
    if (lsm1 == NULL || lsm2 == NULL || lsm1->rows != lsm2->rows||
        lsm1->cols != lsm2->cols) {return ;}

    // 创建结果矩阵
    // LB_Sparse_Matrix* result = (LB_Sparse_Matrix*)malloc(sizeof(LB_Sparse_Matrix));
    // lb_sparse_matrix_initn(result, lsm1->rows, lsm2->cols);
    if(lsm1->rows!=result->rows||lsm1->cols!=result->cols)
    {
        lb_sparse_matrix_resize(result, lsm1->rows, lsm2->cols);
    }
    double* temp_value=(double*)malloc(sizeof(double)*lsm2->cols);

    for (int i = 0; i < lsm1->rows; i++)
    {
        memset(temp_value,0,sizeof(double)*lsm1->cols);
        for(int k=0;k<lsm1->rowSizes[i];k++)
        {
            SparseMatrixEntry entry = lsm1->data[i][k];
            temp_value[entry.i]+=entry.v;
        }

        for(int k=0;k<lsm2->rowSizes[i];k++)
        {
            SparseMatrixEntry entry = lsm2->data[i][k];
            temp_value[entry.i]-=entry.v;
        }
        result->rowSizes[i]=0;

        for(int j=0;j<lsm1->cols;j++)
        {
            if(fabs(temp_value[j])>1e-10)
            {
                result->rowSizes[i]++;
            }
        }
        SAFE_FREE(result->data[i]);
        result->data[i]=(SparseMatrixEntry*)malloc(
            sizeof(SparseMatrixEntry) *result->rowSizes[i]);
        int temp_sum=0;
        for(int j=0;j<lsm2->cols;j++)
        {
            if(fabs(temp_value[j])>1e-10)
            {
                result->data[i][temp_sum].i=j;
                result->data[i][temp_sum].v=temp_value[j];
                temp_sum++;
            }
        }
    }
    free(temp_value);
    // return result;
}


LB_Sparse_Matrix* lb_sparse_matrix_transpose(LB_Sparse_Matrix* lsm)
{
    if (lsm == NULL)
    {return NULL;}

    // 创建新的稀疏矩阵用于存储转置结果
    LB_Sparse_Matrix* transposed = (LB_Sparse_Matrix*)malloc(sizeof(LB_Sparse_Matrix));
    lb_sparse_matrix_initn(transposed, lsm->cols, lsm->rows);

    // 计算每一列的元素数量
    for (int i = 0; i < lsm->rows; i++)
    {
        for (int j = 0; j < lsm->rowSizes[i]; j++)
        {
            SparseMatrixEntry entry = lsm->data[i][j];
            transposed->rowSizes[entry.i]++;
        }
    }

    // 为每一列的数据分配空间
    for (int i = 0; i < transposed->rows; i++)
    {
        lb_sparse_matrix_setrowsize(transposed, i, transposed->rowSizes[i]);
        transposed->rowSizes[i] = 0; // 重置计数器
    }

    // 填充转置矩阵的数据
    for (int i = 0; i < lsm->rows; i++)
    {
        for (int j = 0; j < lsm->rowSizes[i]; j++)
        {
            SparseMatrixEntry entry = lsm->data[i][j];
            int col = entry.i;
            int index = transposed->rowSizes[col]++;
            transposed->data[col][index].i = i;
            transposed->data[col][index].v = entry.v;
        }
    }

    return transposed;
}
// template<class T>
// template<class T2>
// int SparseMatrix<T>::SolveSymmetric(const SparseMatrix<T>& M,const Vector<T2>& b,const int& iters,Vector<T2>& solution,const T2 eps,const int& reset){
//     Vector<T2> d,r,Md;
//     T2 alpha,beta,rDotR;
//     Md.Resize(b.Dimensions());
//     if(reset){
//         solution.Resize(b.Dimensions());
//         solution.SetZero();
//     }
//     d=r=b-M.Multiply(solution);
//     rDotR=r.Dot(r);
//     if(b.Dot(b)<=eps){
//         solution.SetZero();
//         return 0;
//     }

//     int i;
//     for(i=0;i<iters;i++){
//         T2 temp;
//         M.Multiply(d,Md);
//         temp=d.Dot(Md);
//         if(temp<=eps){break;}
//         alpha=rDotR/temp;
//         r.SubtractScaled(Md,alpha);
//         temp=r.Dot(r);
//         if(temp/b.Dot(b)<=eps){break;}
//         beta=temp/rDotR;
//         solution.AddScaled(d,alpha);
//         if(beta<=eps){break;}
//         rDotR=temp;
//         Vector<T2>::Add(d,beta,r,d);
//     }
//     return i;
// }


void lb_sparse_matrix_mult_array(LB_Sparse_Matrix* lsm, double* vec,double* result) {
    if (lsm == NULL || vec == NULL) {
        return ; // 输入无效
    }

    // 分配结果向量空间
    // double* result = (double*)malloc(sizeof(double) * lsm->rows);
    // for (int i = 0; i < lsm->rows; i++) {
    //     result[i] = 0; // 初始化结果向量的每个元素为 0
    // }
    memset(result,0,sizeof(double)*lsm->rows);
    // 遍历稀疏矩阵的每一行
    for (int i = 0; i < lsm->rows; i++)
    {
        // 遍历当前行的非零元素
        for (int j = 0; j < lsm->rowSizes[i]; j++)
        {
            SparseMatrixEntry entry = lsm->data[i][j];
            // 矩阵元素与向量元素的乘法
            result[i] += entry.v * vec[entry.i];
        }
    }

    // return result;
}

void lb_sparse_matrix_printself(LB_Sparse_Matrix *lsm) {
    if (lsm == NULL) {
        printf("Sparse Matrix is NULL.\n");
        return;
    }

    printf("Sparse Matrix:\n");
    printf("Rows: %d, Cols: %d\n", lsm->rows, lsm->cols);

    // 假设 n_nz 是非零元素的总数
    int n_nz = 0;
    for (int i = 0; i < lsm->rows; i++) {
        n_nz += lsm->rowSizes[i];
    }
    printf("Non-Zero Elements: %d\n", n_nz);

    // 遍历每一行
    for (int i = 0; i < lsm->rows; i++) {
        SparseMatrixEntry *current_entry = lsm->data[i];
        for (int j = 0; j < lsm->rowSizes[i]; j++) {
            printf("(%d, %d ,%lf) ", i, current_entry[j].i, current_entry[j].v);
        }
        printf("\n");
    }
}
static inline void conjugate_gradient_compute_residual(LB_Sparse_Matrix* A,
    double* b,double* x,double* ri,double* temp)
{
    lb_sparse_matrix_mult_array(A,x,temp);

    for(int i=0;i<A->rows;i++)
    {
        ri[i]= b[i]-temp[i];
    }

}

static inline double conjugate_gradient_compute_alpha(LB_Sparse_Matrix* A,
    double* r,double* p,double * temp)
{
    double re=inner_product(r,r,A->rows);

    lb_sparse_matrix_mult_array(A,p,temp);
    // //printf("end mult aary\n");
    re/= inner_product(p, temp,A->rows);

    // //printf("alpha:%lf\n",re);
    // free(temp);
    return re;
}

void compute_conjugate_gradient(LB_Sparse_Matrix*lsm,double* b,double * xi ,int iter,double eps)
{

    // double* xi=(double*)malloc(sizeof(double)*lsm->rows);
    double* ri=(double*)malloc(sizeof(double)*lsm->rows);
    double* pi=(double*)malloc(sizeof(double)*lsm->rows);
    double* temp_array=(double*)malloc(sizeof(double)*lsm->rows);
    memset(xi,0,sizeof(double)*lsm->rows);
    memset(ri,0,sizeof(double)*lsm->rows);
    memset(pi,0,sizeof(double)*lsm->rows);
    memset(temp_array,0,sizeof(double)*lsm->rows);

    // conjugate_gradient_compute_residual(lsm,b,xi,ri);
    memmove(pi,ri,sizeof(double)*lsm->rows);
    double threshold=inner_product(ri,ri,lsm->rows);

    for(int i=0;i<iter;i++)
    {
        // printf("begin alpha\n");
        double alpha= conjugate_gradient_compute_alpha(lsm,ri,pi,temp_array);
        for(int j=0;j<lsm->rows;j++)
        {
            xi[j]+=alpha*pi[j];
        }
        // printf("begin residual\n");
        conjugate_gradient_compute_residual(lsm,b,xi,ri,temp_array);
        //printf("end residula\n");
        double temp=inner_product(ri,ri,lsm->rows);
        double beta=temp/threshold;
        threshold=temp;
        // printf("begin pi\n");

        for(int j=0;j<lsm->rows;j++)
        {
            pi[j]=beta*pi[j]+ri[j];
        }


        if(threshold<eps){break;}
    }


    SAFE_FREE(ri);SAFE_FREE(pi);SAFE_FREE(temp_array);
    // return xi;
}
void compute_lu_decompose(LB_Sparse_Matrix*lsm,
    LB_Sparse_Matrix* P,LB_Sparse_Matrix* L,LB_Sparse_Matrix*U)
{
    lb_sparse_matrix_resize(P,lsm->rows,lsm->cols);
    lb_sparse_matrix_resize(L,lsm->rows,lsm->cols);
    lb_sparse_matrix_resize(U,lsm->rows,lsm->cols);

}
